TRANSPARENT BOUNDARY CONDITIONS BASED ON THE POLE 
CONDITION FOR TIME-DEPENDENT, TWO-DIMENSIONAL 

PROBLEMS 
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Abstract. The pole condition approach for deriving transparent boundary conditions is ex- 
tended to the time-dependent, two-dimensional case. Non-physical modes of the solution are iden- 
tified by the position of poles of the solution's spatial Laplace transform in the complex plane. 
By requiring the Laplace transform to be analytic on some problem dependent complex half-plane, 
these modes can be suppressed. The resulting algorithm computes a finite number of coefficients 
of a series expansion of the Laplace transform, thereby providing an approximation to the exact 
boundary condition. The resulting error decays super-algebraically with the number of coefficients, 
so relatively few additional degrees of freedom are sufficient to reduce the error to the level of the 
discretization error in the interior of the computational domain. The approach shows good results 
for the Schrodinger and the drift-diffusion equation but, in contrast to the one-dimensional case, 
exhibits instabilities for the wave and Klein-Gordon equation. Numerical examples are shown that 
demonstrate the good performance in the former and the instabilities in the latter case. 
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1. Introduction. Transparent boundary conditions (TBCs) are required when- 
ever a problem is posed on a domain that has to be truncated in order to become 
numerically treatable, either because it is unbounded or too large to compute solu- 
tions in a reasonable amount of time. Usually, TBCs have to avoid reflections at the 
artificial boundary, although more complex situations can arise, for example if inho- 
mogeneities are present in the truncated part. Exact TBCs are typically non-local in 
time and space and suitable approximations have to be derived in order to be able to 
efficiently compute numerical solutions to the truncated problem. The study of this 
type of boundary conditions started in the 1970s, see the paper of E. L. Lindman [15] 
and references given there. In their seminal paper [5] Engquist and Majda devised a 
general strategy for the derivation of approximate TBCs. Comprehensive overviews 
of the subject can be found, for example, in [21 [71 151 1§1 125] . 

The pole condition approach for the derivation of TBCs was introduced in a first 
version in [2IT1[22] for time-dependent Schrodinger-type equations, later in [T21 IT51 fT§] 
for time-harmonic scattering problems. It was further explored in [6j [TTJ [16j [2T] . An 
alternative formulation of the pole condition is presented in [17J . which provides a 
noticeably simplified implementation and is also used in the present paper. A com- 
parison of different techniques to derive TBCs for Schrodinger's equation can be found 
in [2J, finding the pole condition to be one of the most efficient. In [TH], the pole condi- 
tion approach is adopted for a larger class of time-dependent problems, showing good 
performance for different types of partial differential equations (PDEs) ranging from 
Schrodinger's equation, the heat equation to wave and Klein-Gordon equation. How- 
ever, the experiments involved only one-dimensional or two-dimensional wave-guide 



'Institute of Computational Science, USI Lugano, CH-6904 Lugano, Switzerland. 
E-mail: daniel .ruprechtOusi . ch. Supported by the Swiss HP2C initiative. 

^Mathematisches Institut, Heinrich-Heine-Universitat, D-40255 Diisseldorf, Germany. 
E-mail: schaedle@am.uni-duesseldorf .de 

*ZIB Berlin, D-14195 Berlin, Germany. 
E-mail: frank.schmidtazib.de. Supported by the DFG Research Center Matheon "Mathematics 
for key technologies" in Berlin. 



2 



D. Ruprecht, A. Schadle, F. Schmidt 



geometries. The present paper extends this approach to the fully two-dimensional 
case and investigates its performance through numerical experiments. While the very 
good performance of the pole condition is confirmed in the two-dimensional case for 
Schrodinger's equation and the drift-diffusion equation, instabilities are found for the 
wave equation. 

As the infinite element method, see [3], the pole condition does not truncate the 
exterior domain at some finite length. Nevertheless, the finite number of expansion 
coefficients of the Laplace transform also results in some form of truncation and the 
pole condition realizes a radiation boundary condition at the boundary of the in- 
terior domain and does not aim at providing a meaningful solution in the exterior. 
In some special cases, see [25], the pole condition is closely related to the perfectly 
matched layer approach introduced in [3], but as it does not require complex coordi- 
nate stretching, the pole condition provides a more general framework. Note that in 
contrast to other approaches to TBC involving Laplace transforms, for example [T], 
the pole condition applies the Laplace transform in space and not in time. 

The class of problems considered are, as in |18j . initial value problems for linear 
PDEs of the form 

p(d t )u(t, x) = c 2 Au(t, x) d ■ Vu(t, x) - k 2 u(t, x) for xeR 2 . t> 0. (1.1) 

Included here are the Klein-Gordon equation for p{d t ) = d u and d = (0,0) T , the 
drift-diffusion equation for p(dt) — dt and k — 0, the heat equation for p(di) = dt and 
d = (0, 0) T , k — and finally Schrodinger's equation for p(dt) = idt and d = (0, 0) T , 



k = 0. Equation ( |1.1[ ) is to be solved on a finite computational domain fl C M 2 with 
some boundary condition B(u) = on dfl, such that on the domain Q the solution 
of the initial boundary value problem approximates the solution of the unrestricted 
initial value problem. 

If the support of the initial value u(0,x) is a subset of and the exterior domain 
is homogeneous, in the linear case the boundary condition has to suppress all modes 
traveling from the exterior IR 2 \f2 into the computational domain. Section [2] illustrates 
the main concept of the pole condition by means of a simple one-dimensional example. 
Section [3] introduces the details of the discretization employed in the two-dimensional 
case and Section [4] shows several numerical examples. 

2. Pole condition. This section provides a brief sketch of the key idea of the 
pole condition. Denote the Laplace transform of some function / along some (spatial) 
coordinate r by 



C(f)(s) = / exp(-sr)/(r) dr. (2.1) 
Jo 

The pole condition exploits the identity 

c 1 

cxp(ar) H> , (2.2) 

s — a 

that is a mode with phase a in physical space corresponds to a pole of the Laplace 
transform located at a. The poles of the Laplace transform of the solution are decom- 
posed into poles corresponding to incoming and outgoing modes or, more generally, 
into poles corresponding to physical and non-physical modes. If the locations of the 
poles in the complex plane corresponding to these modes can be separated by a line, 
one can decompose the complex plane into a half-plane Ci n containing all incoming 
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Table 2.1 

Region Ci n for different equations as derived in fT8$. 



Equation 


Parameters in (1.1) 


C in 


Schrodinger equation 
Drift-diffusion equation 

Wave equation 
Klein-Gordon equation 


p(d t ) =id t , d = 0,k = 

P (d t ) = d u k = o 
p(d t ) = d tt ,d = o,k = o 

p{d t ) = d tt , d = 


{z€C: Re(» > -Im(z)} 
{z£C: Re(z) > 0} 
{z E C : Im < 0} 
{z E C : Im < 0} 



modes and a half-plane C ou t containing all outgoing modes. Note that these half- 
planes depend on the equation at hand: Table [2~l] quotes the regions corresponding 
to the equations mentioned above as derived in |18j . For a given Ci n , the pole condi- 
tion is then defined as follows: 

Definition 2.1. Let u(t,r) be a function depending on time t and some spatial 
coordinate r. Denote its Fourier transform in t by u and the dual variable to t by 
lu. Then u satisfies the pole condition, ifU(cu,s) :— £({t(w, -))(s) has an analytic 
extension to Ci n for every u>. 

To illustrate this concept, consider the one-dimensional wave equation on a semi- 
infinite interval 

d tt u{x,t) = d xx u(x,t), x € Cl = [— l,oo) (2-3) 
and assume that a boundary condition at x = is sought such that the solution of 



(2.3) coincides with the solution on the restricted domain [—1,0]. Here, the r from 



Definition |2.1| is identical to the spatial coordinate x. Inserting an ansatz 

u(x, t) — exp(— itot) exp(ifcx) (2-4) 



into (2.3) yields the dispersion relation 

uj = ±k (2.5) 

and assuming uj > without loss of generality yields solutions of the form 

u(x, t) = ci exp (—iujt) exp (iujx) + ci exp (—iujt) exp (— iuix) , (2-6) 

where the first term corresponds to the positive branch of the dispersion relation and 
is rightward moving while the second term corresponds to the negative branch and is 
leftward moving. Let the non-physical modes in this example be the modes traveling 
leftwards from (0, oo) into the interval [—1, 0]. The pole condition then has to suppress 



the pole corresponding to the second term in (2.6) 



In order to point out the connection between the two modes in (2.6) and their 



corresponding poles, we derive an equation for U = C(u) from (2.3). The Laplace 
transform satisfies the identity 

C(d xx f){s) = s 2 £(f)(s) - sfo - f , (2.7) 

where / and /q denote the Dirichlet and Neumann data at x = 0. Further, as in 
Definition |2.1[ denote by U the function obtained by applying to u Fourier transform 



in time and Laplace transform in space. Using ( 2.7 1 , we obtain from ( 2.3 1 the equation 
—uj 2 U (u>, s) = s 2 U(uj, s) — sup(uj) — u (lu) 
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where uq(lo), Uq(u>) are the Fourier transforms of the time-dependent Dirichlet and 



Neumann data uo(t), u' Q (t) at x = 0. By (2.2 1, the first term with pole at iu) cor 



responds to the physically correct rightward propagating mode with coefficient C\ 



in (2.6), the second term with pole at — iuj to the non-physical leftward traveling 
mode with coefficient c^. In order to exclude the non-physical mode, one could for 
example set Qn ={ztC: Im(z) < 0}, so that the pole condition requires U to be 
analytic at — iu>, thus removing this pole from U and the corresponding mode from 
the solution. Note that in this simple one-dimensional example, the pole condition 



can also be enforced by requiring the numerator of the right term in (2.8) to vanish, 
leading to 

iuuo + u' = <H- d t u Q — u' = 0, (2.9) 

which is the well known transparent boundary condition for the one-dimensional wave 
equation, see for example [S]. 



For more complex problems, an explicit decomposition of U like (2.8) is usually 
not available. However, the Laplace transform can often still be decomposed into 
incoming and outgoing parts in terms of path integrals in the complex plane, see |18j , 
hence still allowing to define a region C; n and a pole condition based transparent 
boundary condition. In particular, this is possible for the different types of equations 
listed in Table O 



3. Discretization. This section presents the employed discretizations. Subsec- 
tion |3 . 1 1 describes the discretization of the exterior domain with special semi-infinite 
elements and how the pole condition is incorporated. The mapping between the 
exterior elements and the corresponding reference element introduces a generalized 
distance coordinate, along which the pole condition is enforced. As the exterior el- 
ements are semi-infinite, integrals arise that have one limit infinite. A mapping is 
introduced, converting these integrals into proper integrals in the Hardy space on the 
complex unit-disc. The discretization of the interior uses standard finite elements and 
is not elaborated further. Subsection |3.2| describes the employed integration schemes, 
including the choices of the parameter of the mappings. 

3.1. Space discretization. The discretization of the exterior domain uses semi- 
infinite trapezoids as proposed in The construction of these exterior meshes is 
discussed in [14] . Possible other choices would be semi- infinite triangles and rectangles 
as in [IT] . In any case the exterior discretization has to be such that there is a uniform 
distance variable. For the sake of simplicity, we assume that the computational domain 
f2 is convex, although a generalization to star-shaped domains should be possible. 

3.1.1. Transformation of exterior elements. We mainly use the notation 



of [T7]- Figure 3.1 sketches the used mesh including the exterior elements and the 
mapping between the semi-infinite trapezoids in the exterior and the corresponding 
reference element. Each of the semi-infinite elements T is the image of the semi- 
infinite reference rectangle [0, 1] x [0, oo) under the bilinear mapping g for a set of 
parameters (h„, h^, a, b). Denoting the Jacobian matrix of g by J and its determinant 
by | J |, the mass and stiffness integrals of the variational formulation, which may by 
found in |17j , are repeated below for the convenience of the reader. Additionally the 
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Fig. 3.1. Employed basic mesh. Elements triangulating the interior domain are marked in 
dark grey while the trapezoids decomposing the exterior are marked in light grey. Finer meshes are 
generated by uniform refinements of the triangles and adding additional rays and exterior elements 
when new nodes on the boundary emerge. 




Fig. 3.2. Transformation mapping the reference semi-infinite rectangle to an exterior trape- 
zoidal element. 



drift term is given. 



d ■ V x u v dx 



[0,l]x[0,oo] 



J- T V^ii- J- T Vrtv\J\d(r],0, 



[0,1] X [0,oo ] 



d J- T V rli uv\J\d(r,,0- 



u v dx 



u v\J\ d(r),£), 



'[0,l]x[0,oo] 

The Jacobian of the bilinear mapping g and its determinant are 



J = R 



"" + ( o a + 6)C - b+ t e + b)V )> \J\ = k(h v + (a + b)0, 



(3.1) 



(3.2) 



where R is a rotation, h v is the width of the trapezoid, is a scaling factor measuring 
the distance to the boundary and a and b are signed distance variables, compare for 
Figure |X2] 

As in [T7] we use a product ansatz u(rj,^) :— u(g(?7, £)) — u^^u^r]) on the 
reference strip, where and u v are functions in £ and rj, respectively. The stiffness 
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term, for example, is given by 

f J- T V ni u-J- T V vt :v\J\d(r),Z) 

./ [0,1] X [0,oo ] 

1 ,00 / q ~ - \ T ( hl + (b-{a+b)r,) 2 b-{a+b)r, \ / a ~ ~ \ ( 3 - 3 ) 

d n u v u £ \ [ h7l +e(a+b) — Si — 1 ( d V V V V i 

{a+b)r) h n +{c 



JO 



h„+£(a+6) h ( Vn u V^ d£ , 

U n d £ U £ J I b-(a+b) v fe, + (a+fe)£ I I V n d£Vt ' ? 



Functions u v will be approximated using standard finite element basis functions 
{4>j{vi)}, j = 1,---,N V and therefor integrals over the radial coordinate 77 can be 
evaluated by quadrature formulas. 

3.1.2. Hardy Space. Infinite integrals over the radial coordinate £ are trans- 
formed to finite integrals in the Hardy space H + (Dq) using the identity 

mmn = -2s ±- [ (McmmcgizWzi (3.4) 



2?r JdD 

where M. denotes the modified Mobius transform 

H~(P S0 ) -> H + {D ) : F h> MF defined by (MF)(z) := F (s ^^\ 

(3-5) 

Here, P So denotes a half-plane in the complex plane, depending on the parameter Sq, 
and Do denotes the complex unit-disc, see Figure [3~3] Further, M. is an isomorphism 
between the Hardy spaces H~(P So ) and H + (Dq). Details can be found in (IB ) 117 ] . 

The parameter sq has to be chosen such that the half-plane P So coincides with 
the half-plane C; n of non-physical poles for the considered problem. As functions in 
the space H~(P Sg ) are analytic on the half-plane P So and A4 is an isomorphism, a 
function MF is analytic on Do if and only if F is analytic on P SQ . 

Hence the pole condition, stating that the Laplace transform C(f) of some func- 
tion / has to be analytic on C; n , is equivalent to the condition that AiC(f) is analytic 
on Dq for the correct choice of the parameter sq. In short, we established the following 
sequence of reformulations of the pole condition 

/ satisfies pole condition F has analytic extension to Ci n 

F G H~ (P Sa ) for correct choice of Sq (3-6) 
<^MF e H + (D ), 

see [TB] for details. 

In order to derive a formulation, which is easy to implement, some more transfor- 
mations are required: Given a function /, its image under C and M. is decomposed 
into 



MC{f){z) = — (/o + (z - l)F(z)) =: — T- ™ j (z) (3.7) 

in order to get a local coupling with the boundary data fo. As the Laplace transform 
maps differentiation to multiplication by Sq(z + 1)/(z— 1), straightforward calculation 
yields 

MC{f)(z) = i(/ + (z + l)F(z)) =: T + ( f ° \ (z). (3.8) 
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Fig. 3.3. Sketch of the Mobius transform M : H~(P ag ) —> H~*~(Dq), mapping functions defined 
on the complex half-plane P so to functions defined on the complex unit-disc Do- The parameter so 
has to be chosen such that P SQ corresponds to the Cj n suitable for the problem at hand. 



For motivating this coupling, note that it follows from general theory on Laplace 
transforms that if F is the Laplace transform of / one has lim s _ J . 00 sF(s) — /(0), 
whenever /(0) exists. For details on the decomposition we again refer to [TBI 117] . 

It remains to take care of terms which contain multiplications by £ and ((a + 
b)£ + hr,)- 1 in p3| . To this end, an additional operator V : H+(D Q ) -> H+(D Q ) 
(multiplication by £) Qis implicitly defined by 

ML{{-)f{-)} = M{- (£/(•))'} = s^V (MCf) . 
Direct calculations yield 



(VF) (z) 



(z-lf 



F'(z) + 



F e H+{D 



(3.9) 



Assembling the discrete system involves basically the assembly of discrete counterparts 
of the operators 7+, T- and V. 

3.1.3. Choosing Basis Functions. Up to now, all transformations were on a 
continuous level and no approximations were made. Because functions that are ana- 
lytic on the complex unit-disc can be expanded in power series, the set of monomials 
{V}^_ constitutes a basis of H + (D ). By using the space spanned by a finite num- 
ber of monomials {^'}^ as test and ansatz space, one obtains finite dimensional 
approximations of 7± 



/ 1 



T 



1 / 



-1 
1 



-1 



1 / 



(3.10) 



and V 



Pi 



-1 
1 



V 



(3.11) 



(N S - 1) -(2JV { + 1) J 



The Hardy space monomials can also be transformed back to give a representation of 
the corresponding ansatz and test functions in physical space, see [16] . To define the 



1 This operator is denoted by D in |17l . 
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local stiffness matrix, set for the ^-integrals 



L 



(-i) ._ 



-2Tl _ (Vo I+(a + b)P N( ) 'V, 



(o) 



'2 Tjv ei _IV 4) +, £^2i 



(o) 



-2 Tm i Ta 



.12 — - ' *A t .— J -i-21 •— " - L N i ,+- L N i , 

(!) ._ of, tT t r(°) 



(3.12) 



^5,22 : ~ _2 S 7 W 5 , + 7 A r 5 ,+ ' ^,22 : ~ ~ 2 ( fl + b ) T N (: + P N ( T N ^ + . 

Here the superscript counts the leading order in sq and the subscripts correspond to 



the position in the matrix in (3.3). For the ^-integrals set 

Lr,,U ■= ( I <t>'M (h + 

£77,12 := 
£77,21 := 
£77,22 : = 



((a + b)r) - b)' 



<l>i(v) t Mv)) 

tl£ '1,3=1 

1 b-(a + b)r) Al , n". 
W ^ (Pjiv) 



<i>'M 



(3.13) 



1 



tl£ / i j = l 



For equations with second order temporal derivatives, the parameter s is chosen to 
be frequency dependent in Fourier space, to be precise Sq = iuj, translating back to 
dt in physical space. To avoid the inversion in ^ in (3.12), additional unknowns 
are introduced such that the local stiffness matrices are given by 



£ 



(o) _ 



£77,22 ® £^°22 + £>j,i2 <8 £^12 + £77,21 ® £^21 — 2£ ni n (gi 



-(0) 



(0) 



2I®T N( _ -21 ® (a + b)M Ni 



^77,22 ® £^22 







-2/ (g) /i^7 



(3.14) 



Similarly, local mass matrices corresponding to the mass integral in (3.1 ) are given by 

(3.15) 



M, 



A-i) 



M { ( 1] ®M V 











M, 



(-2) 



Af ? ( 2) ®M n 











where 



M, 



(-1) 



M, 



(-2) 



-2h(h v T N(i _T}f ( - 
-2hz(a + b)Tl _P N T N 



and M v := 



1 N A 

4>i(v)(l>j(v) 



i,3 = l 



(3.16) 



For the drift term set 



n(°) — _?r T 

^£,1 - _ Z1 Nt:,- 



D 



:=-2(a + &)^ ! _P A , 5 T 7V£!+ 



(3.17) 



£,3 • 
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and 



D rhl := h v d 2 [ I 4>i(r])(j)j(r]) 



D V ,3 



= «2 



1 *J=1 



(3.18) 



Mv)(d2(b - (a + b)ri) + dxh^'Ar}) 



i,j=l 



where (di,d,2) T = d = Rd is the rotated d vector. The local drift matrices are then 
given by 



^loc ■ — 



Dfl ® D v<1 



(3.19) 



In the computational domain fi standard local finite element matrices M^l , LrJ 
without the sp-parameter are obtained. By assembling the local matrices to global 



matrices, a spatial semi-discretization of (1.1 1 is obtained 



j(d t ) (V°) + 1 m(-» + ^ m( ~ 2) ) «(*) = ( i(0) + s o^ (1) ) «(*)+ 

^D(o) + il)^ 1 ^ u(t) - k 2 |j|f( ) + 1m + ^2 M(_2) ) «(*) 



(3.20) 



where u(t) is the time-dependent vector of degrees-of-freedom, including the coeffi- 
cients of the monomial basis functions of the subset of H + (Dq) providing the boundary 
condition. 

3.2. Time discretization. All conducted simulations rely on the method-of- 
lines approach: The PDE at hand is first discretized in space, as described in sub- 
section 3.1 leading to the ODE (3.20) for the coefficients. This equation is then 



integrated in time using different time-stepping schemes indicated below. 



3.2.1. Schrodinger's equation. For Schrodinger's equation, (3.20) is solved by 



the second order accurate, A-stable trapezoidal/mid-point rule. Denoting by u n the 
approximation to u{nh) at t = nh for some time-step size h, the discretization reads 



d t Mu(t) = -ic 2 Lu(t) <-» M- 



,n+i 



-i<fL- 



(3.21) 



where the mass and stiffness matrix are given 

and L = + SqL^ 1 '. As in the one-dimensional case, non-physical solutions corre- 
spond to poles in the first quadrant, hence so is chosen in the third quadrant and set 
to 



•so 



-1 



(3.22) 



in order to exclude poles in the region Cj n indicated in Table 2.1 
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3.2.2. Drift-Diffusion equation. In the examples for the drift-diffusion equa- 



tion, (3.201 is integrated with the A-stable Radau II A method with three stages of 
order five, see Sec. IV. 5]. The parameter Sq is chosen to be real and negative, such 
that the poles in the positive half-plane are excluded, corresponding to non-physical 
exponentially increasing solutions. We set 

s = -5, (3.23) 



hence poles with positive real part, see Table pO| are excluded. The chosen value of 
so produces good results, but some optimization is probably still possible. However, 
the sensitivity of the results to the specific value is rather low, as long as the correct 
half-plane is excluded. 

3.2.3. Wave equation. For the wave equation, poles in the lower complex half- 



plane have to be excluded, see 2.1 For the wave equation it is p(dt) — dtt, correspond- 
ing to p(uj) = —to 2 in frequency space. As in the one-dimensional case, we choose the 
parameter so to be frequency dependent, setting 

s = (3.24) 

Transforming back to physical space yields 

d tt M^u - dtM^u + M^u = L^u - d t L^u. (3.25) 

Discretization is done again with the implicit trapezoidal rule, resulting in 

M (Q) U ~ lu + u M ( - 1] — I M ( - 2) = 
h 2 2h 4 

r(0) U ' l+1 + 2u " + r(l) + 1 



2h 

(3.26) 

4. Numerical results. The computational domain in all simulations is a square 
[—4, 4] x [—4, 4] in the two-dimensional plane with slightly smoothed corners, see 
Figure [3~T| Sketched in light gray are the trapezoidal elements spanned by the rays 
in the exterior domain while the darker triangles correspond to the triangulation of 
the interior domain. In order to obtain higher resolutions, the shown mesh is refined 
using up to five uniform refinement steps. As the original mesh is very coarse, no 
errors are reported for simulations on the unrefined grid, because at least lower order 
finite elements do not produce reasonable solutions there. 



4.1. Schrodinger's equation. For p(dt) = idt, c = k = and d = (0,0) 



T 



equation (1.1) yields Schrodinger's equation. For this case, exact solutions of the 
form 

i f-i(x 2 +y 2 )-a(x + y)-2a 2 t\ 

u a (x,y,t;a) = exp (4.1) 

y ' y ' ' ; 4t + i \ U + i J K ' 

with a parameter a are available. We use a superposition of two such solutions, that 



u(x, y, t) = u a (x, y,t;a = 1.4) + u a (x, y, t; a = -2), (4.2) 

and employ u{x, y, 0) as initial value. The discretization in space employs finite el- 
ements of orders one to four in the interior. Simulations are run until T = 10 with 
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Fig. 4.1. Verification of the spatial order of convergence for the Schrddinger equation. Shown 
is the maximum of the relative l^-errors over all outputs versus the refinement level of the mesh 
for finite elements of order one to four. The employed time-step is At = 1/6000 and Ne = 20 
coefficients along each ray are used. As a guide to the eye, lines with slopes one to four have been 
added. 



time-step 


loglO(error) 


Conv. Rate 


1/800 


-3.8 




1/1600 


-4.4 


2.0 


1/3200 


-5.0 


2.0 


1/6000 


-5.5 


2.0 



Table 4.1 

Maximum relative li-error over all generated outputs for the Schrddinger equation, depending 
on the time-step size. The simulation used finite elements of order four, a five times refined mesh 
and iVj = 20 coefficients along each exterior ray. 



time-steps At = 1/800, 1/1600, 1/3200, 1/6000 on meshes refined up to five times and 
for values of (coefficients per ray) between N% = 1 and N% = 20. Output is gener- 
ated at two hundred points in time, distributed equally over the time interval [0, 10]. 



Figure |4.1| shows the maximum relative ^-error over all generated outputs versus 
the refinement levels of the mesh for finite elements of order one to four. All elements 
converge with the expected rate or better until the error saturates at about 3 x 10~ 6 in 
the case of the higher order finite elements. At this point, the error from the temporal 
discretization starts dominating, compare for Table 4.1 and increasing the accuracy 
of the spatial discretization yields no more improvement unless the accuracy of the 
time-discretization is also increased. 

Table [4~T] shows the maximum relative /2-error versus the length of the time-step, 
confirming the convergence rate of two expected from the employed second order 
accurate trapezoidal rule. 
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Fig. 4.2. Relative I2- error depending on the number of coefficients per ray in the exterior 
domain for Schrddinger's equation. Maximum error over all generated outputs for finite elements 
of order one to four on a three times refined mesh. 



Figures 4.2 and 4.3 show how the error decays with increasing N%. The former 
shows the maximum relative ^-error for simulations with finite elements of order one 
to four, a three times refined mesh and a time-step of A = 1/6000. In all cases, the 
error decays super-algebraically with N% until it saturates at the level of the respective 
spatial or temporal discretization error. Figure |4.3| shows the error at four different 
points in time for the fourth order elements. The error decays super-algebraically at 
all four points in time, even at later times where most of the solution has left the 
domain. 

Figure 4.4 shows the relative ^-error over time for different values of N%. In all 



cases, the error increases as the wave packets hit the boundary of the computational 
domain and subsequently decays to a level determined by the number of coefficients 
per ray N%. The error decays faster for larger values of N^, but for values of N$ = 10 
or larger, the levels at which the error saturates and in particular the error at the end 
of the simulation is identical. 

4.2. Drift-diffusion equation. Setting p(dt) — dt and k = in yields the 
drift-diffusion equation. Note that the heat equation is included here as the special 
case d = (0, 0) T . An analytic solution is given by 



u(x, y,t) = - exp 



1 

4te2 



[x - dttf + (y- d 2 tf 



(4.3) 



Set 



d = (di,d 2 ) = (1.5,1.5), and c = 0.5 



(4.4) 



and start the integration at to = 0.2 with initial value u(x, y, to). This yields a 
Gaussian function with a peak initially close to the origin which is subsequently 
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Fig. 4.3. Relative l^-error at four fixed points in time depending on the number of coefficients 
per ray for the simulation with finite elements of order four. In all cases, a time-step of A = 
1/6000 has been used. 



time-step 


loglO(error) 


Conv. Rate 


1/10 


-3.6 




1/20 


-5.0 


4.6 


1/30 


-5.9 


4.9 


1/40 


-6.5 


4.9 


1/60 


-7.4 


4.9 


1/80 


-8.0 


5.0 


1/160 


-9.2 


4.0 



Table 4.2 

Maximum relative l^-error over all generated outputs depending on the time-step size for finite 
elements of order six, a four times refined mesh and Nc = 51 coefficients per ray. 



advected to the upper right corner of the square while being spread out by diffusion. 
Integration in time is done by the fifth order Radau IIA(5) scheme. The simulations 
are run until T — 5 with finite elements of order one to six, on meshes refined up to 
four times, values of N% between 1 and 51 and time-steps ranging from A = 1/10 to 
At = 1/160. 

Figure 4.5 shows the maximum relative ^-error over all generated outputs versus 
the refinement level of the mesh. The number of coefficients per ray is N$ — 51 and 
the time-step is At = 1/160. As a guide to the eye, lines with slopes from one to six 
are added. All elements converge with the expected rate or better, confirming again 
that the pole condition does not compromise the order of convergence of the spatial 
discretization. 

Table |4.2| shows the maximum relative ^-error versus the length of the employed 
time-step for finite elements of order six, a four times refined mesh and N% — 51. 
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Fig. 4.4. Relative l^-error over time for the Schrodinger equation for six different values of 
Ne . The simulation used finite elements of order four, a five times refined mesh and a time-step 
At = 1/6000. 



From At — 1/10 to At = 1/20, a slightly reduced convergence rate is observed, 
probably because the time-step size is still in the pre-asymptotic regime. The reduced 
convergence rate in the last refinement is because the error approaches the spatial 
discretization error, see Figure |4~5| Beside that, the expected fifth order convergence 
is observed, demonstrating that the pole condition can not only preserve the accuracy 
of high order finite elements but also of high order integration schemes. 



Figures 4.6 and 4.7 show the relative ^-error versus the number of coefficients 
ATj . The former shows the maximum error over all outputs for finite elements of order 
one to six, a three times refined mesh and a time-step At = 1/160. In contrast to 
Schrodinger's equation, for small values of there is only a minor decrease of the 
error. Also, below N^, the error is not decreasing monotonically with the number of 
coefficients. After N$ = 10, rapid super-algebraically decrease of the error is again 
observed until the error saturates at a level determined by the accuracy of the spatial 
discretization. Note that N% = 30 coefficients per ray are sufficient here to reduce 
the boundary condition error to the level of the spatial discretization error in all 
cases. Figure |4.7| shows the relative error at four different points in time. Again the 
error generally decays super-algebraically with the number of coefficients, but now 
the decay rates are noticeably lower at later points in time. Note that at t = 0.43, the 
Gauss peak has not yet reached the boundary, so that the boundary condition has no 
visible effect on the error at this time. 



Figure 4.8 shows the relative ^-error over time for different values of JVjt, finite 
elements of order six and a four times refined mesh. As for Schrodinger's equation, 
the error starts increasing at some point in time, but the increase starts later as the 
number of coefficients increases. On the other hand, for the simulations with N% = 15 
or less coefficients, the errors at the end of the simulation are about the same and 



Transparent boundary conditions 



15 




Fig. 4.5. Verification of the spatial order of 
grated until T = 5. Shown is the maximum of the 
number of coefficients per ray is = 51 and the 



convergence for the drift- diffusion equation inte- 
relative l^-errors over all generated outputs. The 
time-step is At = 1/160. 



only for larger values of N$ a significantly reduced error is observed at the end of 
the simulation. Together with Figure |4~6{ this illustrates that for the drift-diffusion 
example, the error is not monotonically decreasing with N% for small values of N$ and 
a certain minimum number of coefficients per ray is required before the onset of the 
super-algebraic decay. 

4.3. Wave and Klein-Gordon equation. For p(dt) = du and d = (0, 0) T , 



( 1.1 ) becomes the Klein-Gordon equation, containing the wave equation as the special 
case k = 0. While the pole condition could successfully provide TBCs for both 
equations in the one-dimensional case as well as in a two-dimensional wave-guide 
problem, stability problems arise in the fully two-dimensional case, rendering the pole 
condition in the here presented form inapplicable to both equations for finite elements 
of order two or higher. Resolving these issues is planned for future research. 
Below, the instability is documented briefly. Use an initial distribution 

u{x,y)=exp(-2x 2 - 2y 2 ) , (4.5) 

finite elements of order one to four and up to four refinement steps for the mesh. 
Integrate in time using implicit trapezoidal rule with a time-step At = 1/1280 until 
T = 100. 

Figure [4~9| shows the maximum ^-error over all generated outputs versus the size 
of the elements of the employed mesh for finite elements of order one to four. In 
order to obtain a readable plot, the error is capped at 10 1 and values of 10 correspond 
to unstable runs. While first order elements are stable on all five meshes, showing 
again better than expected convergence, for higher order elements and fine meshes, 
the method becomes unstable. Second order elements are unstable only on the finest 
mesh while third and fourth order elements already become unstable after three or 
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Fig. 4.6. Relative l2-error depending on the number of coefficients per ray in the exterior 
domain for the drift-diffusion equation. Maximum error over all generated outputs for finite elements 
of order one to six on a three times refined mesh. 



two refinement steps. Note that on the coarser meshes where the method is stable, 
the expected or better decay rates of the error are observed. Similar behavior is found 
for the Klein-Gordon equation, but not documented here. 



Figure 4.10 shows the energy of the discrete solution over time for different values 
of N% on the finest mesh for finite elements of order two. The simulations are stable 
for Nt = 4. They are also stable until about t = 20 for AT £ = 6, 10, 12, 15, but an 
exponential instability occurs after this point in time. The instability also occurs when 
the simulation is run on a circular domain. Note that the employed integration scheme 
is A-stable, so that the instability on a finer mesh is not arising from a violation of 
some CFL-type stability limit. 

5. Conclusions. The pole condition approach to transparent boundary condi- 
tions, derived in |18j for the time-dependent, one-dimensional case, is extended to 
time-dependent two-dimensional problems. The pole condition identifies in- and out- 
going modes by associating them with poles of the spatial Laplace transform in the 
complex plane. The complex plane is then divided into two half-planes, C; n and 
C ou t, containing the poles corresponding to incoming and outgoing modes respec- 
tively. To suppress modes traveling from the exterior into the computational domain, 
the Laplace transform is required to be analytic in Cj n . In order to obtain a numer- 
ically implementable formulation, Ci„ is mapped to the unit circle by a conformal 
Mobius transformation. The Laplace transform is then extended in a power series on 
the unit circle with the coefficients of the expansion being connected to the interior 
degrees of freedom on the boundary. Truncating the series after a finite number of 
terms yields an approximate and implementable TBC. 

Numerical examples are presented in order to investigate the performance of the 
pole condition approach: As in the ID-case, the considered generic PDE contains 
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Fig. 4.7. Relative l^-error depending on at four fixed points in time for the simulation with 
finite elements of order six. In all cases, a time-step At = 1/160 has been used. 



different well-known equations for specific choices of parameters. Excellent results are 
obtained for Schrodinger's equation and the drift-diffusion equation: The presented 
numerical experiments demonstrate that the convergence order of finite elements up to 
order six is retained and also that the convergence order of the temporal discretization 
is not affected if sufficiently many coefficients are used for the boundary condition. 
Further it is shown that the error introduced by the approximate boundary condition 
decays super-algebraically as the number of coefficients in the expansion of the Laplace 
transform increases. For the drift-diffusion equation, a small minimal number of 
coefficients was found to be required to reach the regime of super-algebraic error 
decay. 

Unfortunately, in contrast to the one-dimensional case, the approach exhibits 
instabilities for the two-dimensional wave and Klein-Gordon equation if using finite 
elements of order two or higher. Hence in the present form the pole condition is of 
limited use for these second order hyperbolic equations. A further investigation of the 
instability and hopefully a remedy will be subject of future research. 
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Fig. 4.9. Spatial order of convergence for the wave equation integrated until T = 100. Shown 
is the maximum l^-error over all generated outputs. The error is capped at 10 1 , so depicted error 
values of 10 correspond to unstable runs. 
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Fig. 4.10. Energy over time for the wave equation for different values of N$. The refinement 
level of the mesh is two and the order of the used finite elements is four. 



